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We briefly review analytic approximations of thermodynamic functions of fully ionized nonideal electron- 
ion plasmas, applicable in a wide range of plasma parameters, including the domains of nondegenerate and 
degenerate, nonrelativistic and relativistic electrons, weakly and strongly coupled Coulomb liquids, classical 
and quantum Coulomb crystals. We present improvements to previously published approximations. Our code 
for calculation of thermodynamic functions based on the reviewed approximations is made publicly available. 
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1 Introduction 

In a previous work (T]21, we performed hypernetted chain (HNC) calculations and proposed analytic formulae 
for the equation of state (EOS) of electron-ion plasmas (EIP). An alternative analytic approximation for the EOS 
of EIP was proposed in |[3]|4|. A comparison (e.g., UMIi) shows that the formulae in have a higher accuracy, in 
particular for the thermodynamic contributions of ion-ion and ion-electron correlations in the regime of moderate 
Coulomb coupling. Recently ||5j|6l, we studied classical ion mixtures and proposed a correction to the linear 
mixing rule. In this paper we review the analytic expressions for all contributions to thermodynamic functions and 
introduce some practical modifications to the previously published formulae. The reviewed analytic description 
of the thermodynamic functions of Coulomb plasmas is realized in a publicly available computer code. 

Let Tie and rij be the electron and ion number densities, A and Z the ion mass and charge numbers, respectively. 
The electric neutrality implies rie = Zni. In this paper we neglect positrons (they can be described using the 
same formulae as the electrons; see, e.g., Ref. Q) and free neutrons (see, e.g., Ref. ISj), and consider mainly 
plasmas containing a single type of ions (for the extension to multicomponent mixtures, see ||5]|6l). 

The state of a free electron gas is determined by the electron number density and temperature T. Instead 
of He it is convenient to introduce the dimensionless density parameter = a^/ao, where qq is the Bohr radius 
and Ge = (iTrne)"^/"^. The parameter can be easily evaluated from the relations = 1.1723 7124 where 
71,24 = 7ie/10^* cm^'^, or rg = {po/pY^^, where po = 2.6752 (A/Z) g cm^^. The analogous density parameter 
for the ion one-component plasma (OCP) is Rs = a{m{{Z ef' / fp' = 1822.89 VgAZ'^/^, where toi is the ion mass 
and Oi = (|7rni)~^/'^ = OeZ^/'^ is the ion sphere radius. 

At stellar densities it is convenient to use, instead of r^, the relativity parameter Q a;r = pp/meC = 
1.00884 {peZ/A^^^^ 0.014005 rj^ where pp = hi^Tr^n^y^^ is the electron Fermi momentum and pe = 
p/10^ g cm"''. The Fermi kinetic energy is ep = c (mgc)^ + pp — nieC^, and the Fermi temperature equals 
Tp = ep/fes = Tr (7r - 1), where Tj. = nieC^ /ks = 5.93 x 10^ K, 7^ = ^1 + x'^, and ks is the Boltzmann 
constant. If <C 1, then Tp sa 1.163 x 10^ r j^ K. The effects of special relativity are controlled by in 
degenerate plasmas (T <C Tp) and by r = T /T^ in nondegenerate plasmas (T S> Tp). 

The ions are nonrelativistic in most applications. The strength of the Coulomb interaction of ions is char- 
acterized by the Coulomb coupling parameter V — (Ze)^ /aiksT = TeZ^^^, where Fg = e? /aeksT sa 
22.747r6~^ [peZ/A)^'^ and Tg = T/IO^ K. 
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Thermal de Broglie wavelengths of the ions and electrons are usually defined as Aj = (2'Kh? /riiikBT) " and 

1/2 

Ae = {2nh^ /m^kBT^ ■ The quantum effects on ion motion are important either at Ai > Ci or at T ^ Tp, 
where Tp = (hujp/kB) ~ 7.832 x 10^ (Z/A) ^/pQ K is the plasma temperature determined by the ion plasma 

1 /2 

frequency cjp = (47re^ TiiZ^/mj) . The corresponding dimensionless parameter is 77 = Tp/T. 

Assuming commutativity of the kinetic and potential energy operators and separation of the traces of the 
electronic and ionic parts of the Hamiltonian, the total Helmholtz free energy F can be conveniently written as 
F = i^/^'' + f//' + Fee + Fa + F,e, where F.j-^^ and f//^ denote the ideal free energy of ions and electrons, 
and the last three terms represent an excess free energy arising from the electron-electron, ion-ion, and ion- 
electron interactions, respectively. This decomposition induces analogous decompositions of pressure P, internal 
energy U, entropy S, the heat capacity Cy, and the pressure derivatives xt = {d\nP/d\nT)v and Xp ~ 
— (9 In P/d In y )t • Other second-order functions can be expressed through these ones by Maxwell relations. 



2 Ideal part of the free energy 

The free energy of a gas of Ni = riiV nonrelativistic classical ions is F^^ = NiksT [ln(?iiAf/A/) - l] , where 
AI is the spin multiplicity. In the OCP, it can be written in terms of the dimensionless plasma parameters as 
^id^ = N-.ksT [3\nrj - l.SlnT - 0.51n(6/7r) - InM - 1] . 

The free energy of the electron gas is given by = ^e^e ^ Pj^^ V, where fie is the electron chemical 
potential. The pressure and the number density are functions of fie and T: 



P 



(e) _ 8 ksT 



T 1 4 

^3/2(Xe,T) + -/5/2(Xe,T)J , Tie = [/l/2 (Xe , t) + T/3/2 (Xe , t)] , (1) 



where Xe = tJ-e/ksT (here, we do not include the rest energy meC in fie) and 

r x'^{l + rx/2y/\ 
Jo exp(a; - Xe) + I 

is a Fermi-Dirac integral. An analytic approximation for fig (n^) has been derived in IT]. 

In Ref. Ill we calculated Iu{Xe: t) using analytical approximations Q. These approximations are piecewise: 
below, within, and above the interval 0.6 < Xe < 14. Their typical fractional accuracy is a few xlO^'', the 
maximum error and discontinuities at the boundaries reach ~ 0.2%. For the first and second derivatives, the 
errors and discontinuities lie within 1.5%. At Xe > 14 we use the Sommerfeld expansion (e.g., ifTOl ) 

/. (xe, r) « 4"^ (A) + y r^I^^) (fi) + ^r^lf^) (fi) + . . . , (3) 
where we have defined fi ~ XeT — fJ-e/rUeC^, 

Jo Jo Vl + x'' dfJ- 

xi^\e) = + e and xo = \/ /i(2 + jl). In particular, 

l[%{fi) = [xo7o-ln(xo + 7o)]/2, = xl/3-l[%{fi), lfj>^{fL) = xi^7o/4- 2x^/3 + 1.25lj°^(/i), 



where 70 = \J\ + x^ = 1 + fi (note that, if /i = e, where e = ep/meC^, then Xo = Xr and 70 = 7r). 

At small fi, accuracy can be lost because of numerical cancellations of close terms of opposite signs in Eq. Q 
and in the respective partial derivatives. In this case we use the expansion 

^ ('2777 - 1 V T™+1 r°° dr 

/.(x,r)=/r(x)+E(-l)'" 4Ji';„, ^"Wi(x), C(X)=/ -^r^, (5) 
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where (2m — 1)!! = nfeLi(2fc — 1) should be replaced by 1 for m = 0. At large the nonrelativistic Fermi 
integrals /""^(x) calculated with the use of the Sommerfeld expansion. Sufficiently smooth and accurate 
overall approximations for functions I,y{xe,T) with = 5> |> and f are provided by switching to Eq. ^ at 
Xe > 14 and fi < 0.1, while retaining the terms up to m ~ 3. 

The chemical potential at fixed rie is obtained with fractional accuracy ^ T'^/Tp by using Eqs. ([T), (O, setting 
li"'\fi) ~ li"'\e) + (// — e), and dropping the terms that contain the product {jl — e) t^. Then at 

where Fq'^'' = (Pi.V'/87r^)[a;r fl + 23;^) "/,■ — Infxr +'7rl] — N^m^c^ is the zero-temperature limit ifTll (without the 
rest energy NeirieC^), and = mgC^ {mf.c/Kf' = 1.4218 x 10^^ dyn cm~^ is the relativistic unit of pressure. 
Accordingly, p!^f' w P^"'^ + AP, where Po*"^ = (Pi /Stt^) [x^ (f - l) 7,. + \-a{x, + 7^)] , and AP = 

(Pr/18) r^x,. (2 + a;2)/7,. In this case, x^'^ « P,x5/97r27,p/;\ Xt^ « 2AP/p/;\ « ^^^'siV, 7rr/x?. 
In order to minimize numerical jumps at the transition between the fit at Xe < 14 and the Sommerfeld expansion 
at Xe > 14, we multiply the expressions for AP by empirical correction factor (1 + Ae/e)~^, and those for Xt^ 
and dy'^ by [1 + (4 - 2a;,./7r) l^e/e]-\ 

At Xi- < 10^^ we replace Pq'^^ and Pq*^^ by their nonrelativistic limits, F^'^^ /V ~ Pi-x^/IOtt^ and Pq'^^ — 
PrX^/lSTT^ cx TT-e^''. In the opposite case of x-^ ^ 1, one has Pq'^' = Pra;J?^/127r^ cx rie^^. 



3 Electron exchange and correlation 

Electron exchange-correlation effects were studied by many authors. For the reasons explained in fTJ, we adopt 
the fit to Pee presented in Ref. 112]. It is valid at any densities and temperatures, provided that ^ 1. 

In Ref. II2I we implemented an interpolation between the nonrelativistic fit |[T2l and approximation fSlIll valid 
for strongly degenerate relativistic electrons. However, later on we found that such an interpolation may cause 
an unphysical behavior of the heat capacity in a certain density-temperature domain. Therefore we have reverted 
to the original formula |[T2l . taking into account that in applications, as long as the electrons are relativistic, their 
exchange and correlation contributions are orders of magnitude smaller than the other contributions to the EOS 
ofEIP (see, e.g., ilini). 



4 One-component plasma 

4. 1 Coulomb liquid 

For the reduced free energy of the ion-ion interaction fa = Fa/NiksT in the liquid OCP at any values of F we 
use the analytic approximation derived in Ref. In order to extend its applicability range from T ^ Tp to 
T ~ Tp, we add the lowest-order quantum corrections to the Helmholtz free energy d: /f ^ = ?7V24. The 
next-order correction oc h"^ has been obtained in I15J . These corrections have limited applicability, because as 
soon as jy becomes large, the Wigner expansion diverges and the plasma forms a quantum liquid, whose free 
energy is not known in an analytic form. 

A classical OCP freezes at temperature Tm corresponding to F = 175, but in real plasmas Tm is affected by 
electron polarization EIISl and quantum effects (e.g., iflTlffQl ). Hence the liquid does not freeze, regardless of 
T, at densities larger than the critical one. The critical density values in the OCP correspond to Rs ~ 140-160 
for bosons and Rs ~ 90- 1 10 for fermions. In astrophysical applications, the appearance of quantum liquid can 
be important for hydrogen and helium plasmas (see, e.g., Ref. IfTSlI for discussion). 



4.2 Coulomb crystal 

The reduced free energy of the Coulomb crystal is /lat = PiatZ-^Vi^'sP = C'oF+1.5ui?7+/th + ./ah- Here, the first 
term is the Madelung energy (Co ~ —0.9), the second represents zero-point ion vibrations energy {ui « 0.5), 
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Fig. 1 Anharmonic contribution to the reduced free energy 
of an ion lattice as a function of the quantum parameter 
77 — Tp/T for two values of the Coulomb coupling param- 
eter, r = 175 (upper curves) and 1000 (lower curves). A 
comparison of different approximations (see text): Wigner 
expansion (long-dashed lines); an approximation from 1251 
(lOI, dotted lines); the same approximation with the coeffi- 
cients adjusted to Ref. 1211 (short-dashed lines); and present 
interpolation ^ (solid lines). The points with errorbars 
show simulation results from 1251 for F = 1000. 



Fig. 2 Harmonic and anharmonic lattice contributions to 
the reduced heat capacity at F = 175 and 1000. The har- 
monic lattice contribution according to Ref. |20| (dot-dashed 
lines) is compared to the model II 181 (long-short-dash lines) 
and to the anharmonic correction in different approximations 
(see text): a derivative of the lOI fit 1251 (dotted lines: orig- 
inal; short-dashed lines: improved; see text) and a corre- 
sponding derivative of Eq. ^ (solid lines). 



/th is the thermal correction in the harmonic approximation, and /ah is the anharmonic correction. We use the 
most accurate values of Co, ui, and analytic approximations to /th for bcc and fee Coulomb OCP lattices, which 
have been obtained in ll20l . 

For the classical anharmonic corrections, 11 fitting expressions were given in |23]- We have chosen one of 
them: /^°^(r) = ai/T + a2/2T'^ + ag/SF^ with ai = 10.9, 02 = 247, and as = 1.765 x lO'^, because this 
choice is most consistent with the perturbation theory 11711221 . 

In applications one needs a continuous extension for the free energy to ?/ 7^ 0. With the leading quantum 
anharmonic correction at small 1] one has IfTSl 

/ah « /^(r) - (0.0018/r + o.085/r2) r,\ (?) 

At r — > 0, the quantum anharmonic corrections were studied in Ill7ll23ll24l . where an expansion in powers of 

— 1/2 

Rg was obtained, assuming Rs ^ 1. The leading term of this expansion can be written in the form /ah.T^o = 
—birf'/V. In the literature one finds different estimates for 61; we use hi = 0.12 as an approximation consistent 
with 12311241 . Free and internal energies of finite-temperature quantum crystals were calculated using quantum 
Monte Carlo methods in ||19ll25ll26l . lyetomi et al. in Ref. Il25l (hereafter lOI) proposed an analytic expression 
for the quantum anharmonic corrections. However, differences between the numerical results in |fT9ll25ll26ll 
are comparable to the differences between the numerical results and the harmonic approximation. Therefore, 
finite-temperature anharmonic corrections cannot be accurately determined from the listed results. 

In order to reproduce the zero-temperature and classical limits, we multiply fj^^ by exponential suppression 

2 

factor e~'^^^ and add /ah.T-s-o- Then we have 

/ah = /r(r)e-^^'''-6ir;Vr. (8) 
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According to Eq. the Taylor expansion coefficient at 77^ /F equals zero. In order to reproduce this property, 
we set ci — bi/ai « 0.0112. In Fig.[T] the resulting approximation for the free energy as a function of 77 (solid 
lines) is compared to the semiclassical expression [Eq. Q; long-dashed lines] and to the elaborated lOI fitting 
formula: the dotted lines correspond to the original coefficients of this fit, chosen by lOI so as to reproduce the 
results of Ref. Il22l at 77 — !• 0, while the short-dashed curves show the same fit, but with the coefficients adjusted 
to more recent results 1211 at 77 — )■ 0. The curves of the two latter types nearly coincide at F = 1000, but differ 
near the melting point F = 175. 

Interpolation (|8j (unlike, for example, Pade approximations ll27l ) does not produce an unphysical behavior of 
thermodynamic functions. In particular, anharmonic corrections to the heat capacity and entropy do not exceed 
the reference (harmonic-lattice) values of these quantities at any 77. In Fig.|2] we show the ion contributions to the 
reduced heat capacity of the ion lattice cv.i = Cv.i/NikB, calculated through derivatives of Eq. (O (solid lines), 
and compare it to the contributions calculated through derivatives of the fitting formula in i25'| (short-dashed 
and dotted lines). In the same figure we plot the harmonic-crystal contribution to the reduced heat capacity: the 
short-dash-long-dash line corresponds to the model of Ref. ifTSl (which we adopted in ||2|), while the dot-dashed 
line corresponds to the most accurate formula 1201 . 

Interpolation Eq. (O should be replaced by a more accurate formula in the future when accurate finite- 
temperature anharmonic quantum corrections become available. 



5 Electron polarization 

5 . 1 Coulomb liquid 

Electron polarization in Coulomb liquid was studied by perturbation 11311281 and HNC |[l]|2]|29]|30l techniques. 
The results for Fi^ are fitted by analytic expressions in |2j|. This fit is accurate within several percents, and it 
exactly recovers the Debye-Hiickel limit for EIP at F — and the Thomas-Fermi result 19] at F 1 and Z ^ 1. 



5.2 Coulomb crystal 

Calculation of thermodynamic functions for a Coulomb crystal with allowance for the electron polarization is a 
complex problem. For classical ions, the simplest screening model consists in replacing the Coulomb potential by 
the Yukawa potential. For instance, there were molecular dynamics simulations of classical Yukawa systems ||T6] 
and path-integral Monte Carlo simulations of a quantum Yukawa crystal at Rs — 200 ISTl . However, the Yukawa 
interaction reflects only the small-wavenumber asymptote of the electron dielectric function. A rigorous treatment 
would consist in calculating the dynamical matrix and solving a corresponding dispersion relation for the phonon 
spectrum. The first-order perturbation approximation for the dynamical matrix of a classical Coulomb solid with 
the polarization corrections was developed in Ref. Il32l . The phonon spectrum in such a quantum crystal has been 
calculated only in the harmonic approximation ll33l . which has a restricted applicability: for example, it cannot 
reproduce the known classical (T ^ Tp) limit of the anharmonic ie contribution to the heat capacity. 

A semiclassical perturbation approach was used in ||2l- The results were fitted by the analytic expression 

U = -fooix,) F [1 + A{x,.) iQ{v)/m ■ (9) 



In the classical Coulomb crystal, Q = 1, foo{x) = 61 \/l + and A{x) = (63 + a3x'^)/{l + b/^x^). 

Parameters s and 61-64 depend on Z and are chosen so as to fit the perturbational results; 03 is a constant. All 
the parameters weakly depend on the lattice type; they are explicit in [2j for the bcc and fee lattices. 

The factor (5(77) in Eq. (|9]i is designed to reproduce the suppression of the dependence of Fi^ on T at T <C Tp. 
For the classical solid, Q(0) — 1. In the quantum limit (r/ 00) Q{ri) ~ q-q, so that the ratio Q/T in Eq. (|9]) 
becomes independent of T. The proportionality coefficient q was found numerically in [121; it equals 0.205 for 
the bcc lattice. 

The form of Q{ri) suggested in ||2l assumed a too slow decrease of the heat capacity {Cv,ie oc 77"^ oc T) at 
?7 — >^ 00, which signals the violation of the employed semiclassical perturbation theory in the strong quantum 
limit, related to the approximate form of the ion structure factor used in the calculations, as discussed in IJIEI. In 
order to fix this problem, we have changed the form of Q{ri) to 

-1/2 

(10) 



1 /2 

Q(77) = [in (1 + e**")')] [in (e - (e - 2)e-('"')') 
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This form of Q{i]) has the correct limits at 77 and 77 ^ 00 and is compatible with the numerical results IJ)- 
In addition, it eliminates the problematic ie contributions to the heat capacity and entropy at 7] ^ 1. It can 
be improved in the future when the polarization corrections for the quantum Coulomb crystal are accurately 
evaluated. 



6 Conclusions 

We have reviewed analytic approximations for the EOS of fully ionized electron-ion plasmas and described 
recent improvements to the previously published approximations, taking into account nonideality due to ion-ion, 
electron-electron, and electron-ion interactions. The presented formulae are applicable in a wide range of plasma 
parameters, including the domains of nondegenerate and degenerate, nonrelativistic and relativistic electrons, 
weakly and strongly coupled Coulomb liquids, classical and quantum Coulomb crystals. 

For brevity we have considered plasmas composed of electrons and one type of ions. Extension to the case 
where several types of ions are present is provided by [5 6 1 . 

We have made the Fortran code that realizes the analytical approximations for the free energy and its deriva- 
tives, described in this paper, freely available in the Internet^. 
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